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Bubble  oscillations  of  large  amplitude 

Joseph  B.  Keller  and  Michael  Miksis 

Departments  of  Mathematics  and  Mechanical  Engineering,  Stanford  University,  Stanford ,  California  94305 
(Received  7  December  1979;  accepted  for  publication  21  April  1980) 

A  new  equation  is  derived  for  large  amplitude  forced  radial  oscillations  of  a  bubble  in  an  incident  sound  field. 
It  includes  the  effects  of  acoustic  radiation,  as  in  Keller  and  Kolodner's  equation,  and  the  effects  of  viscosity 
and  surface  tension,  as  in  the  modified  Rayleigh  equation  due  to  Plesset,  Noltingk  and  Neppiras,  and 
Poritslty.  The  free  and  forced  periodic  solutions  are  computed  numerically.  For  large  bubbles,  such  as 
underwater  explosion  bubbles,  the  free  oscillations  agree  with  those  obtained  by  Keller  and  Kolodner.  For 
small  bubbles,  such  as  cavitation  bubbles,  with  small  or  intermediate  forcing  amplitudes,  the  results  agree 
with  those  calculated  by  Lauterbom  from  the  modified  Rayleigh  equation  of  Plesset  et  al.  For  large  forcing 
amplitudes  that  equation  yielded  unsatisfactory  results  whereas  the  new  equation  yields  quite  satisfactory 
ones. 

PACS  numbers:  43.30.Jx,  43.30.Gv,  43.30.Lz,  43.25.  Yw 


INTRODUCTION 

A  gas  bubble  in  a  liquid  performs  forced  radial  os¬ 
cillations  when  a  sound  wave  impinges  upon  it.  Large 
amplitudes  result  when  the  acoustic  frequency  is  at  or 
near  the  bubble’s  natural  frequency,  or  certain  rational 
multiples  of  it.  Since  these  large  oscillations  can  be 
important  in  cavitation,  we  have  reexamined  them  by 
using  a  new  equation  which  takes  account  of  the  acoustic 
radiation  from  the  bubble.  The  previous  detailed  and 
careful  calculations  by  Lauterborn'  did  not  take  account 
of  this  radiation,  and  they  sometimes  yielded  unrea¬ 
sonably  large  amplitudes  which  did  not  lie  on  smooth 
curves,  or  they  failed  to  converge  to  a  periodic  solu¬ 
tion.  Our  method  does  not  have  these  defects. 

Oscillations  of  large  bubbles  were  originally  analyzed 
by  Rayleigh  assuming  that  the  surrounding  liquid  is  in¬ 
compressible  and  inviscid,  that  the  bubble  remains 
spherical,  and  that  surface  tension  is  negligible.  Ples¬ 
set,2  Noltingk  and  Neppiras,3  and  Poritsky'  modified 
this  equation  to  include  the  effects  of  viscosity,  surface 
tension,  and  an  incident  sound  wave,  and  it  was  this 
modified  equation  which  Lauterborn1  solved.  A  differ¬ 
ent  modification  was  made  by  Keller  and  Kolodner,5 
who  included  the  effects  of  acoustic  radiation  by  treat¬ 
ing  the  surrounding  liquid  as  slightly  compressible. 

The  same  method  was  used  by  Epstein  and  Keller0  to 
derive  equations  for  one-  and  two-dimensional  bubbles, 
for  which  there  are  no  analogues  of  the  Rayleigh  equa¬ 
tion.  We  shall  combine  these  modifications  to  derive 
a  new  equation  for  the  bubble  radius.  It  includes  the  ef¬ 
fects  of  acoustic  radiation,  viscosity,  surface  tension 
and  an  incident  sound  wave. 

To  solve  this  new  equation  we  shall  first  study  its 
free  oscillations.  As  is  to  be  expected,  the  trajectories 
in  the  phase  plane  are  similar  to  those  in  the  absence 
of  viscosity  and  surface  tension  unless  the  viscosity  is 
very  large.  The  oscillations  of  a  particular  underwater 
explosion  bubble  are  found  to  be  the  same  as  those  cal¬ 
culated  by  Keller  and  Kolodner,5  which  agree  well  with 
experiment. 

Then  we  solve  the  equation  numerically  for  one  of 
the  three  bubbles  studied  by  Lauterborn,1  with  various 
forcing  amplitudes.  Our  results  are  in  very  close 


agreement  with  his  for  small  and  intermediate  forcing 
amplitudes,  but  not  for  his  largest  one.  In  that  case 
his  results  suffer  from  the  defects  mentioned  above, 
while  ours  are  satisfactory  for  this  and  even  for  larger 
forcing  amplitudes.  We  have  also  examined  the  small 
amplitude  oscillations  by  the  method  of  averaging,  fol¬ 
lowing  the  procedure  used  by  Prosperetti7  on  the  equa¬ 
tion  of  Plesset  et  al.,  but  we  shall  not  present  that  an¬ 
alysis  here. 

I.  FORMULATION 

We  wish  to  determine  a(l).  the  radius  at  time  /  of  a 
spherical  bubble  of  gas  and  vapor  in  a  fluid  of  infinite 
extent.  We  assume  that  the  pressure  />„(«)  within  the 
bubble  is  uniform  and  is  the  sum  of  the  constant  vapor 
pressure  pv  and  the  gas  pressure  ka~*r: 

/>,(«)  =  />„  ♦fen-”.  (1.1) 

Here  it  is  a  constant  determined  by  the  quantity  and  type 
of  gas,  and  y  is  the  adiabatic  exponent  of  the  gas.  This 
pressure  pb{a)  exceeds  the  pressure  p(r  .1)  in  the  fluid 
at  r-a  by  the  effect  of  the  surface  tension  a  and  the 
normal  component  of  viscous  stress: 

/>,(«) =/>(«,/)+ ^2-^ ^rr(a,t)  (1.2) 

Here  <p(r,t)  is  the  velocity  potential  of  the  fluid  motion, 
assumed  to  be  spherically  symmetric,  and  p  is  the  co¬ 
efficient  of  viscosity  of  the  fluid.  In  addition  the  velo¬ 
city  a,  of  the  bubble  surface  must  equal  the  velocity  of 
the  fluid  at  the  surface: 

a,-  <br(a,t) .  (1.3) 

The  velocity  potential  0,  the  pressure  p.  and  the  den¬ 
sity  p(r,l)  of  the  fluid  must  satisfy  conservation  of 
mass,  the  Navier-Stokes  equation,  and  the  equation  of 
state, which  are 

p,+ <prpr+ pA<t>=  0,  (1.4) 

p(0r,+  0r0rr)  +  pr=  (4(1/3)(A0),  ,  (1.5) 

p=p(p).  (1.6) 

We  seek  a,  p,  />,  and  0  satisfying  (1.1)-(1.6)  given  their 
initial  values. 
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II.  SIMPLIFICATION 

To  solve  these  equations,  we  first  divide  (1.5)  by  p 
and  then  integrate  it  with  respect  to  r  from  >-  to  infinity 
to  obtain  the  modified  Bernoulli  equation 

+ £  p'1  '//>  —  (4/x  3 )J  p-'(±<t>)r,lr  =  0.  (2.1) 

Here  />„  is  the  constant  pressure  at  »•  =  •*>,  and  we  have 
assumed  that  O  tends  to  zero  at  »-=  ■».  Next  we  dilrer- 
entiate  (2.1)  with  respect  to  /  and  write  t>  t-  PpP,  =  >'~Pt. 
where  <'•’  =  /> 0  and  c  is  the  sound  speed: 

+  <V>rl+  '''P|/>''  +  |p',(Ad>)r|,r/>  =  0 .  (2.2) 

Then  we  use  (1.4)  in  (2.2)  to  eliminate  pt.  divide  by  r~ 
and  obtain 

c"26((  -  ±<t) 

=  ~<''2  C'V’rt +  (V  3)J^  l p',(A<.‘>), p  . 

(2.3) 

To  simplify  (2.3)  we  omit  the  right  side,  assuming 
that  it  is  small  compared  to  the  terms  on  the  left  side, 
and  set  c  =  constant.  For  a  nearly  incompressible  fluid 
«•*  is  large  and  pr  is  small,  which  tends  to  justify  these 
assumptions.  As  a  consequence  we  obtain  the  wave 
equation 

0.  (2.4) 

To  simplify  (2.1)  we  set  p=  constant  in  it  and  obtain 

f>(r, /)  =  />. *'!>?!+  (4p  3)A<6  .  (2.5) 

Thus  the  simplification  consists  in  replacing  (1.4)  and 
(1.5)  by  (2.4)  and  (2.5). 

III.  DERIVATION  OF  THE  ORDINARY  DIFFEREN¬ 
TIAL  EQUATION 


=  ^nj+  a,Ma)  1  +  —  -  A(n)\ 

2  *  1  \2  ‘  pa  pa  j 

+  rni,A'(rt)+  2^1  +  ~^^'"^  +  ~y  13.4) 

Equation  (3.4)  is  a  nonlinear  second-order  ordinary 
differential  equation  for  the  bubble  radius  a{l).  If  we 
set  p  =  o  =  A'  =  0  in  it,  we  obtain  the  corresponding  Eq. 

( 1 2)  of  Keller  and  Kolodner  .5  If  instead  we  divide  by  c 
and  letr  become  infinite,  we  obtain  the  Rayleigh  equa¬ 
tion  for  a  bubble  in  an  incompressible  liquid .  as  mod¬ 
ified  by  Plesset2  and  others  to  include  surface  tension 
and  viscosity.  That  is  the  equation  which  was  solved 
by  Lauterborn1  with  2r ■’/<"(/)  replaced  by  p~'P  sinu.v . 

In  order  to  solve  (3.4)  we  must  specify  ^  and  the  initial 
values  of  a  and  at.  Then  (1 .3)  yields/,  (3.3)  yields  cb 
and 2.5  yields  />. 

In  order  to  specify  ^  we  suppose  that  the  bubble  is  in 
an  incident  sound  field  with  velocily  potential  4>(x./), 
where  x  =  0  is  the  center  of  the  bubble.  The  spherically 
symmetric  part  of  4>  is  of  the  form  r’l|A'(/4  >’  r) 

+  //(/-»•  r)|.  If  4>  is  regular  at  x  =  0  then  h  -  and  we 
have 

4>(x./)  =  r‘l|.v(f  +  r  c)-iiU-r  <)|+  4>„(x./) .  (3.5) 

The  unsymmetric  part  4>u  must  vanish  at  x=  0,  so  if  we 
let  x  and  r  tend  to  zero  in  (3.5)  we  gel 

4>(0./)=  2r'V>,(/) .  (3.6) 

This  relation  determines  #'(/)  in  t.  rms  of  the  incident 
potential  <t>  evaluated  at  the  position  of  the  center  of  the 
bubble.  Differentiating  (3.6)  yields 

2*"(/)  =  r4>,(0./).  (3.7) 

This  expression  is  what  is  needed  in  the  last  term  in 
(3.4). 


We  now  use  (2.5)  for  p  in  (1.2)  and  write  that  condi¬ 
tion  in  the  form 

at  r  =  <»(/) .  (3.1) 

Here  A(«)  is  the  pressure  difference  times  p*1: 

A((j)  =  p*,(p>(«)  -/>„ \=  p',[kn"'r+pv-p.\.  (3.2) 

Next  we  write  the  general  solution  of  (2.4),  with/  and  g 
arbitrary  functions,  as 

<t>lrj)=r"'fU  -r  r)*r''nd*r  r).  (3.3) 

Then  we  use  (3.3)  in  (3.1)  to  obtain  an  equation  for  <i(/). 
In  doing  so  we  also  use  (3.3)  in  (1.3)  to  eliminate  f.  In 
this  way  we  obtain 


riAtri)  -  aca,  -  oMn. /)  - 


1 

2 


rmj  + 


2a 

n 


As  an  example,  let  us  suppose  that  the  incident  field 
is  a  plane  wave  with  angular  frequency  u>  and  pressure 
amplitude  P.  Then  4>(x,/)=  -(P  pw) cosur(/  -  v  c)  and 
(3.7)  yields 


=  Pep’1  sinic/  .  (3.8) 

We  now  substitute  (3.8)  into  (3.4)  to  obtain 


a“(^-n{n>  ~r)) 


+  a,M< i)  - 


4  pa,  2a 
pa  pa 


+  imtA'(n)+  ^1  4  sin 4  . 


(3.9) 


This  is  the  equation  we  shall  solve.  In  the  Appendix  it 
is  rewritten  in  dimensionless  variables. 


IV.  ANALYSIS  OF  THE  EQUATION 


(3.3') 

Now  we  take  the  time  derivative  of  (3.3')  and  use  (3.1) 
to  eliminate  .  In  this  way  we  obtain 


When  I*=  0,  (3.9)  reduces  to  the  autonomous  equation 
for  the  free  oscillations  of  a  bubble.  This  equation  can 
be  analyzed  in  the  phase  plane  with  coordinates  a  and  a, . 
It  has  one  unstable  critical  point,  a  saddle,  at  n=0. 
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(i,  =  c  (3}  -  1)  and  a  second  critical  point  at  the  equili¬ 
brium  configuration  a  =  ne,  <i,=  0.  where  ae  is  a  root  of 

2(7  (Wt-a|(!c)=0,  (4.1) 

Linearization  shows  that  this  point  is  a  stable  spiral  if 
u)o>0  and  A  stable  node  if  where  w*  is  given 


5,  using  the  appropriate  parameter  values.  The  results 
are  indistinguishable  from  those  given  in  that  figure, 
which  were  computed  with  a=  p  =  0.  This  was  to  be  ex¬ 
pected  since  the  dimensionless  forms  of  o  and  p  are  so 
small  in  that  case. 

The  periodic  solution  of  the  linearized  form  of  (3.9) 
withP*  0  is 

fl(()-ne  =  ~Iln[.Aexp(iiu(<  +  r'1<7t)]  ■  (4.4) 

where 


Usually  a’J>0  and  the  equilibrium  point  is  a  spiral. 
However  if  the  fluid  is  extremely  viscous  it  will  be  a 
stable  node.  When  u’0  is  real  the  linearized  solution 
near  equilibrium  is,  with  A  and  9  constants. 


a{l)  -ne=A  exp 


cos(u.'0/  +  e ) . 


14.3) 


This  describes  a  damped  oscillation  around  equilibrium. 


For  the  nonlinear  equation  (3.9)  with  P=0.  when  a.’* 

>  0  the  trajectories  in  the  phase  plane  are  similar  to 
those  for  ci=  p  =0,  which  are  shown  in  Ref.  5.  We  have 
also  solved  (3.9)  numerically  with  P  =  0  for  the  case  of 
the  underwater  explosion  bubble  treated  in  Fig.  7  of  Ref. 


As  c  -  A  tends  to  A„  given  by 


(4.5) 


.  /  ,  .  .,  ,  2a  4ptcu\'‘ 

=  (  -a.ur-  - - )  . 

V  P«*  1 


We  can  make  A„  agree  with  A  by  replacing  p  in  (4.6)  by 
a  certain  complex  effective  viscosity.  Alternatively  we 
can  make  [A.  |  agree  with  |A|  by  replacing  p  in  (4.6) 
by  the  real  quantity  /aeI[  defined  by 


FIG.  1.  Frequency  response  curves 
for  a  bubble  in  water  with  equili¬ 
brium  radius  ae  =  10  pm  in  incident 
sound  waves  of  amplitudes  P=  0.4, 
0.7,  0.8,  and  0.9  bar.  The  ordinate 
is  (amax~ae)/ae  and  the  abscissa 
is  The  curves  are  computed 

from  (3.49)  and  u>o  is  given  by  Eq. 
(4.2). 
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a«s«  .  , 


p 

FIG.  3.  Amplitude  of  the  subharmonic  response  as  a  function  of  the  amplitude  of  the  incident  sound  wave  for  w/u>t» 2  and  uj/u>g  >3. 


-pa,^(n,)  +  a,A'{a,)j  ]+  G(c'J)  (4.7) 

An  approximate  form  of  (4.7)  appropriate  for  high  fre¬ 
quencies,  p.„  =  p+ p<tjwj/4c,  was  derived  by  Devin,* 
Chapman  and  Pies  set*  and  Prosperetti.10 

V.  NUMERICAL  SOLUTIONS 

In  solving  (3.9)  numerically  we  shall  follow  the  pro¬ 
cedure  used  by  Lauterborn1  on  the  corresponding  equa¬ 
tion  for  an  incompressible  fluid,  in  order  to  show  that 
the  differences  in  results  are  due  to  differences  in  the 
equations,  and  not  due  to  differences  in  methods  of 
solution.  Thus  we  start  with  a  bubble  at  rest  at  its 
equilibrium  radius,  i.e.  <r,(0)  =  0,<i(0)  =  n,.  Then  for 
given  Pand  u>,  we  follow  the  solution  until  it  becomes 
periodic,  and  plot  (rtmu -n,)/n,  vs  u>.  u>„.  The  constants 
are  chosen  to  represent  a  bubble  in  water  at  20°  C  with 
a  static  pressure  A'l  bar  and  a  polytropic  exponent 
y  »  1.33.  Thus  p*  0.998  g/cm3,  a*  72.5  dyn  cm,  p 
*0.01  g  cm  s,  p„*0.0233  bar  and  c *  1.484  x  10s  cm/s. 
We  shall  study  Lauterborn ’s  largest  bubble,  for  which 


a,=  10'3  cm,  since  acoustic  radiation  is  most  important 
for  it. 

Figure  1  shows  our  frequency  response  diagram  for 
four  values  of  the  incident  pressure  amplitude:  P  =  0.4, 
0.7,  0.8,  and  0.9  bar.  It  corresponds  to  Lauterborn 's 
Fig.  3.  For  P*  0.4  bar  our  curve  and  his  appear  to  be 
identical,  while  for  P *  0.7  bar  they  are  also  identical 
except  for  the  ultraharmonic  resonance  of  order  5/2 
which  occurs  in  his  curve  at  about  u>/u>0  =  0.37 ,  but  is 
absent  in  ours.  To  see  how  this  difference  arises,  we 
shall  consider  the  radius  time  curve  a(l)  for  P*  0.7  bar 
and  ui/w0*  0.37  in  both  the  incompressible  and  com¬ 
pressible  cases.  Figure  2(a)  shows  the  steady  state  or 
periodic  oscillation  in  the  incompressible  case  and 
Fig.  2(b)  shows  some  early  oscillations  in  the  compres¬ 
sible  case.  We  note  that  the  two  curves  are  somewhat 
similar.  However  as  time  goes  on,  in  the  compressible 
case  the  shape  of  the  oscillation  changes  to  the  periodic 
form  shown  in  Fig.  2(c),  which  does  not  contain  the  5  2 
ultraharmonic. 

Our  response  curve  for  P*  0.8  bar  is  very  similar  to 
that  for  P*  0.7  bar,  but  it  is  quite  different  from  Lau¬ 
terborn ’s  curve.  His  incompressible  result  and  our 
compressible  one  are  similar  for  u>/u>(l>0.6,  but  not  for 
u>/u>0<  0.6.  The  incompressible  result  shows  many 
peaks  corresponding  to  harmonic  and  ultraharmonic 
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resonances,  at  some  frequencies  the  amplitudes  were 
too  large  to  fit  in  the  figure,  and  at  some  frequencies 
no  periodic  solution  was  obtained.  The  compressible 
result  always  yielded  periodic  solutions  with  amplitudes 
which  fit  into  the  figure,  and  some  harmonic  and  ultra¬ 
harmonic  resonances. 

For  I’=  0.9  bar  our  result  is  again  similar  to  that  for 
!‘=  0.8  bar.  but  with  more  resonances.  We  computed  it 
only  for  a.'  0.5  to  save  computer  time,  since  the  os¬ 

cillation  went  through  many  cycles  before  becoming 
periodic. 

In  Fig.  1  we  also  see  a  subharmonic  resonance  at  j.1 
j.'„=2.  This  resonance  did  not  occur  for  small  values 
of  P.  but  instead  there  is  a  threshold  value  which  f‘  must 
exceed  in  order  for  it  lo  appear.  This  threshold  phi  - 
nonienon  is  shown  in  Fig.  3.  For  a.-  jl\,  =  2  we  find  that 
the  threshold  value  of  /'lies  between  0.20  and  0.22  Dar 
while  Lauterborn  found  it  to  be  between  0.1  and  0.15 
bar.  Similarly  for  u.'0=3  the  compressible  threshold 
lies  between  1.9  and  1.94  bar  while  in  the  incompres¬ 
sible  case  Lauterborn  found  it  to  be  between  1.3  and 
1.5  bar.  Thus  in  both  cases  the  inclusion  of  compres¬ 
sibility  increases  the  threshold,  as  one  might  expect. 
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APPENDIX 

In  terms  of  the  equilibrium  radius  <>,.  which  satisfies 
4.1).  we  define  dimensionless  variables  as  follows: 

v  =  a  a,,,  r  =  r  nr , 

f>y  'J  "e . 

O  =  f))l/L'.  k-k  />,/!, !r  . 

(Al) 

(f=  2o  (•*<■((>  />„)*' ‘ 

i/  =  4m  <if(p/ >„)l/p  *"  =  (#>  V'  . 


We  now  introduce  these  variables  into  (3.4)  and  then 
omit  the  bars  to  obtain 

v„lr  -i(\,  -r)|=^4+.v([(l  -3y)A>.v'3,-l| 

(Ax:,‘  ex.  a 
~2*~^*x~kx  +l1 


ex.  a  ,  ...  A 

-<(y+T"+x— +v 

*4^MKM£n 


(A2) 


In  the  same  way  we  obtain  from  (4.1) 

A1  =  1  4  a  .  (A3) 

Similarly  we  get  from  (3.9) 

v ,,le  -  v(v(-r)|  =  '^-  +  x,[(l  -3y)kx~sr  -  1 1 
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